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ABSTRACT 

We explore methods of fitting templates to cosmic microwave background (CMB) data, and in 
particular demonstrate the application of the total convolver algorithm as a fast method of performing 
a search over all possible locations and orientations of the template relative to the sky. This analysis 
includes investigation of issues such as chance alignments and foreground residuals. We apply these 
methods to compare Bianchi models of type VII^ to WMAP first year data and confirm the basic 
result of our 2005 paper. 

Subject headings: cosmology: cosmic microwave background - cosmology: observations 



L INTRODUCTION 

The widely accepted model in cosmology, the so-called 
concordance model, posits an isotropic and homogeneous 
universe with small anisotropics generated by primordial 
fluctuations in the infiationary field. These anisotropics 
are present in the cosmic microwave background (CMB), 
which should then be statistically isotropic and Gaussian. 
Many CMB studies therefore examine the CMB from a 
statistical point of view with the intention of testing for 
violations of these properties. Alternative cosmological 
models have not, however, been completely ruled out, 
and there are several anomalies in the Wilkinson Mi- 
crowave Anisotropy Probe ( WMAP) data that indicate 
that such models merit further investigation by alternate 
means. 

We investigate methods for testing any determin- 
istic anisotropic cosmological model. The predicted 
anisotropy template can be compared to the data us- 
ing fitting techniques in both pixel and harmonic space 
to search for correlations. We present a description of 
these methods and apply a fast and efficient algorithm 
for searching the full sky for the best orientation of a tem- 
plate relative to the data. We test these methods with 
both full- and incomplete-sky data sets, and use simula- 
tions to characterize the significance of the results. 

Motivated by the morphology of several detected vi- 
olatio ns of Gaussianity and/or isotropy in the WMAP 
data, llde Olivejra.Costa et al1l2nn3 rF>ik^^;^^ 
iHansen et all l2004bi: iVielva et all I2004D . we test our 
methods using Bianchi type VII^i models and the WMAP 
first-year data. A preliminary analysis was published in 
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iJaffe et all l)2005 () , in which we reported on a surprisingly 
significant detection of a Bianchi model at the 99.7% sig- 
nificance level compared to simulations. Here we present 
an improved search of the model space, confirm the ba- 
sic result, and discuss in detail issues such as foreground 
contamination and chance alignments. 



2. METHODS 



2.1. 



Template Fitting 

Given any anisotropy pattern that contributes to the 
data as an additional component of the observed mi- 
crowave sky (whether topological in origin, as in the case 
of Bianchi models, or foreground), we perform a fit of 
the template to the WMAP data as has been done in 
the past by, e.g., (HsSrski et al. (1996) and Bandav ct aQ 
(1996) for foreground analysis. The best-fit amplitude a 
for a template vector t compared to a data vector d can 
be measured by minimizing 

- (d - at)^M^i (d - at) = d^M.i d, (1) 



where Mgji^ is the covariance matrix including both 
signal and noise for the template-corrected data vector 
d = d — at. Solving for a then becomes 



a = 



t^M 



(2) 



SN' 



To compare multiple template components to a given 
data sets, e.g., different foregrounds, the problem be- 
comes a matrix equation. In the case in which we have 
N different foreground components, we define 



N 

d = d - ^afetfc 

k=l 



and 



(3) 



(4) 



Mgjv^ = (dd"^) = Mg + 

In this case, minimizing d-^Mg-J^^d leads to the following 
set of equations. 



N 



^t^Mg^t.a, =t^Mg^d. 



(5) 
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This is the simple system of hnear equations Ax = b, 
where 

Xk = ak- 

When only one template is present, this reduces to equa- 
tion |j2l above. 

The errors (5a^ are the square root of the diagonal of 
A^^. The matrix A gives information about the cross- 
correlation between the templates themselves. 

Note that the above is equally valid in pixel space or 
harmonic space. In the former, it is very easy to account 
for incomplete sky coverage or to remove, for example 
the Galactic plane, by simply including in the data vec- 
tors only the relevant pixels, and likewise by including 
only the corresponding rows and columns of the covari- 
ance matrix. The noise in pixel space is usually well 
represented by a diagonal matrix representing uncorre- 
lated pixel noise. But the signal covariance matrix in 
pixel space is large and not sparse, which makes har- 
monic space more convenient when this is possible. In 
harmonic space and under the assumption of Gaussian- 
ity, the signal covariance is diagonal, and with the ap- 
proximation of uncorrelated noise that is uniform over 
the sky, the noise covariance can be made to be so as 
well. The difhculty in harmonic space is t he sky cover- 
age. As discussed bv lMortlock et al.l l|20n2D . the coupling 
matrix to cut out the Galactic plane using a cut the size 
of |6| > 20° becomes numerically singular for resolutions 
of ^max > 50. Cuts such as the conservative KpO mask, 
defined by the WMAP team, remove more of the sky 
and, due to their structure, the coupling matrix is more 
difficult to compute. 

We define a method that applies harmonic space fit- 
ting to the full-sky cases using highly processed maps 
discussed in § 13.21 This allows us to increase the co mpu- 
tational efficiency using the algorithm described in § 12.21 
Here, we use a uniform mean noise approximation that 
has a diagonal covariance matrix. We use pixel-space 
fitting for each band separately in the cut-sky analysis 
in which the Galactic plane region is masked out. (At 
the WMAP signal-to-noise ratio level, little would be 
gained by simultaneously fitting the different bands, and 
the memory and CPU requirements to invert the ma- 
trix would become onerous.) Again, we use a diagonal 
noise approximation that this time takes into account 
the observation pattern but not the effects of smoothing. 
Comparisons of fits with fully correct noise to those using 
these approximations show that the results do not vary 
significantly (at most a few percent, or a small fraction 
of the error bar) . All codes have been cross-checked with 
identical inputs to confirm identical outputs. 

Note that the cosmic monopole and dipole are not 
known, and although a best-fit dipole is subtracted from 
the data in the map making process, small residual 
monopole and dipole terms remain in the data. For this 
reason, we cannot include this component in the fit. In 
harmonic space, any monopole and dipole terms can sim- 
ply be excluded or ignored by setting, e.g., Ci = C2 = 
10^/iK^. In pixel space, we fit the monopole and dipole 
simultaneously as independent components. See ii4.3l for 
discussion. 



The above method determines the best-fit amplitude 
for a given template at a fixed orientation relative to the 
sky. For foreground fitting, that is generally all that is re- 
quired, but in the search for an anisotropic cosmological 
component, there are two additional complexities. First, 
we have no a priori guess for where the symmetry axis 
may be pointing and must thus search the entire sky. Sec- 
tion l2.2l describes the algorithm we use to do this quickly 
and efficiently. Second, we may have an infinite number 
of possible templates (e.g., parameterized as described 
in § 13.1(1 among which we want to find the "best" , so in 
addition to determining the best-fit amplitude for each 
template, we need a way to compare how well different 
templates fit the data and to select the most interesting. 
Section discusses how we address this. 

2.2. Total Convolver 

The search for the best orientation of a template com- 
pared to the data requires that we evaluate the statistic 
a described in the previous section at every possible rel- 
ative orientation of the template and data. Working in 
harmonic space allows us to use an algorithm based on 
Fourier transforms to speed up this search significantly. 

In the case of full-sky analysis, the location or orienta- 
tion of the template does not affect the error, i.e., 5a is 
invariant. Then the maximum of a is found at the maxi- 
mum of the numerator in equation |5J above, t-^Mg^d. 
Neglecting for the moment the covariance matrix, the 
quantity to be maximized is simply the convolution of 
the data with the template. We seek the maximum over 
all possible locations and orientations, and this can be 
found effici ently using the total co nvolver algorithm de- 
scribed in Wan delt fc Gorskil l|2001 l). which was originally 
developed for map making using instrument beams. 

This algorithm decomposes the Euler angles into what 
amounts to a scan pattern and then takes advantage of 
the form the convolution takes in harmonic space to sim- 
plify the calculation. The rotation operator -D($2, ©, ^1) 
can be factored into I?((/)jt;, Q)D{(f), 6, oj), where a pre- 
defined scan pattern determines 9-^ and 9, which in the 
case of full-sky coverage are both 7r/2, so that the set of 
angles {(p-p , (f>, u!) covers the full sky at all possible ori- 
entations (see iWandelt fc Gorskil 1(200 ID Figure 1). (In 
only this context of total convolution on the full sky, (j) 
corresponds to the polar angle and (f>j^ to the azimuthal 
angle. Elsewhere in this paper, these are represented by 
the more common 9 and (p.) Defining T{(j)-^, 4>,uj) = t^d 
as the quantity to be maximized, b£m as the spherical har- 
monic coefficients of the template t, and Ofm. as that for 
the d ata d, the convolution is then l(Wandelt fc GorsTil 
l2?)?}Tl eqs. 9 & 8) 

Tmm'm" =^ Q^jmC^Lm' (^E)'^m'm" (6) 
I 

rrt 1 1 \ Sr^ rr im(jiT?+im' 4>+im" uj ,„\ 

T{(j)-^,(p,Uj) = Tram'm"e E , (7) 

where d^^,{9) is the real function such that 

DL^'ih,OAi) = e-;'"^^d^„,(0)e-™''^i. The problem 
has then become simply to calculate Tmm'm" and 
Fourier transform to T((/)g, 0, cj) to find the maximum. 
To take into account the signal and noise covariance, we 
simply use a "whitened" data vector, Mq^d. 
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The total convolver can find the best-fit position with 
an accuracy limited only by the resolution of the inputs. 
The positional accuracy is 7r/£rnax, which for our analy- 
sis is 2°. 8. Note that this is larger than the size of a pixel 
at the usual HEALPix resolution of ^gj^ig = ^niax/2. 

It should also be noted that searching the full sky will 
not return an unbiased estimate for the amplitude. Sim- 
ulations with a known input value for a particular tem- 
plate at a known position will, on average, have slightly 
higher amplitudes returned by the search. If the correct 
template location is simply fit to an ensemble of sim- 
ulations with additional CMB and noise, the returned 
amplitudes will have a Gaussian distribution with the 
correct mean and variance, but the same is not true when 
one is searching for the best location and orientation as 
well. This is because the search is seeking the maximum, 
and the resulting distribution is a form of extreme value 
distribution^, which introduces a small positive bias in 
the results. For realistic situations with CMB and noise 
in addition to the component we arc fitting, the total 
convolver is likely to find a maximum amplitude a small 
distance away from the true position. How different the 
amplitudes and positions are on average depends on the 
particular case in question, since it is a function of how 
dominant the template is compared to the CMB and 
noise, and how much the template structure changes over 
angular distance, etc. This is quantified for the particu- 
lar case in question in M. 21 using simulations. 

This method is approximately 2 orders of magnitude 
faster than performing the fit in harmonic space over 
a grid of individual rotations one at a time. The dis- 
advantage is the storage requirement for the matrix T, 
which increases with the third power of the resolution 
and becomes over 2GB for a HEALPix resolution of 
^side ~ ^"^^ angular resolution of 42'. 

2.3. Best-fit Model and Significance 

As mentioned above, when it is not one unique tem- 
plate for which we are testing but rather a set of possibil- 
ities, we need not only to find the best fit of each to the 
data but also to find the best fit among the possible mod- 
els. Depending on how the model space is parameterized, 
there can be an infinite number of possibilities. Previ- 
ous studies se eking u pper limits on shear and vorticity 
l)Kogut et al.l [T997: Bu nn et al.l ll99d) used two different 
statistics to determine the "best" -fit model. 

Given a model, Kogut et al. define the best-fit position 
and amplitude in terms of F = a/Sa. They used Cos- 
mic Background Explorer (COBE) data, for which no 
full-sky analysis was possible. In the case of incomplete 
sky analysis, the amount of template structure that is 
masked changes the significance of the fit. A large am- 
plitude in which most of the structure is masked by the 
Galactic plane cut is not as interesting as a lower ampli- 
tude fit in which the structure is included. By finding 
the maximum not of a but of F, they attempt to find 
the most significant fit rather than simply the maximum 
amplitude. 

Bunn et al. (1996) use a different statistic to accom- 
plish the same effective selection. They define 771 = 
iXo ~ Xi)/Xo7 where Xi is as in equation and Xo the 

Sec, e.g., " http : //mathworld ■ wolf ram . com/ 1 

ExtremeValueDistr ibut ion . html 



corresponding statistic for the data by itself, uncorrected 
for any anisotropic component. The difference is then 
an indication of how much better the data fit the (sta- 
tistically isotropic CMB) theory after correction for the 
anisotropic model. 

Finding the maximum of F is equivalent to finding the 
maximum of 771 (although Bunn et al. use a different 
statistical method) . So for a given model, either statistic 
can be used to find the best-fit amplitude and position. 
But it becomes more complicated to compare one model 
to another in order to determine which model fits the 
data better. 

The prob lem with the simple approach, used by 
iKogut et all 119971) a s well as in our preliminary anal- 
ysis (• Jaffe et al.ll2005f) , of using F or r/i to find the best 
model is that the distribution of these values for chance 
alignments is not the same for all models. Although they 
are generally quite similar, differences in the tails of the 
distributions mean that a given value of F has a slightly 
different significance for different models. This means 
that finding the maximum of F might have missed other 
models that are significant but in which the tail of the 
distribution does not reach as high in F. In other words, 
the significance of the fit found in our original result is 
not incorrect, but it is possible that such an analysis fails 
to detect another significant model. 

For this more c omplete analysis, we analy ze a set of 
LILC simulations ijEriksen et al.ll2004all2005j) . using the 
above formalism to characterize the distributions of a 
values for a given model. In this analysis, for a given 
model, we compare the a value (equivalently F, since Sa 
does not change for a given model on the full sky) for the 
WMAP data against the ensemble of simulations. We 
can then quantify the significance of a given model fit 
to the data based on the percentage of LILC simulations 
in which the model fits with a lower amplitude. This 
gives clearer indication of which are the most interesting 
models than a simple F or x'^ statistic. Comparison of 
the results using a or 771 in this way show there is little 
difference between the two in terms of how significant a 
given fit to the data is against the simulations. In the 
following analysis, we use the numbers for a only. 

2.4. Visualization: Cross-Correlation Signal Maps 

It is helpful to be able to visualize what parts of the 
sky are driving a particular fit. To do this, we simply 
note that the numerator of equation[51 t-^Mg^d, can be 

rewritten in pixel space as ^p[L~^t]p[L~^d]p, where L 
is the "square root" of the covariance matrix Mg]\j, or 
its lower triangular decomposition found from, for exam- 
ple Cholesky decomposition. A simple visualization is to 
turn this into a map, where each pixel contains the prod- 
uct of [Lg-Jjt] and [Lg-J^d] at that pixel. This map shows 
exactly what regions on the sky drive the fit at a given 
orientation. This is particularly important when certain 
regions of the sky are known to be contaminated; these 
plots show whether or how much those regions affect the 
fit. Examples will be shown in S 15.21 

3. DATA AND SIMULATIONS 

Here we describe the particular class of models we in- 
vestigate and the data sets used in the analysis. 

3.1. Bianchi Models 
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Bianchi type VII/i refers to the class of spatially homo- 
geneous generalizations of Friedmann universes that in- 
clude small vorticity (universal rotation) and shear (dif- 
ferential expansion) components. (Type VIIq includes 
the flat Friedmann-Robertson- Walker model, and VII/i 
include s that with negative spatial curvature as special 
cases.) iBarrow et all 1)19851) solve the geodesic equations 
to derive the induced CMB anisotropy by linearizing the 
anisotropic perturbations about the Friedmann models. 
Their solution does not include any dark energy compo- 
nent, which is a significant shortcoming considering the 
preponderance of evidence that now points to CIa ^ 0.7. 
But we examine them first as a test of our template- 
fitting methods and second because of the intriguing pos- 
sibility that they may explain several anomalies in the 
data. 

Following the prescription in IBarrow et alJ l)1985|) . we 
construct a template for the anisotropy induced by vor- 
ticity (w) and shear (tr). Bianchi type VII/i models are 
parameterized by the current total energy density f2o and 
a parameter x (^Collins fc Hawking 1973^ . 



(8) 



where h is related to the canonical struct ure constants 
and is tha t to whic h the type Vll h refer s (see'Kogu t et all 
[1997; Bun n et al.i ri996: Barrow et all [l985). This pa- 
rameter can be understood as the ratio of the scale on 
which the basis vectors change orientation to the Hub- 
ble radius (present values). The resultin g temperature 
aniso tropy pattern is then described by l)B arrow et al.l 
IT9851 eq. 4.11) 



±[B{e 



R) 



Ai0R)]cos{^R)}, 



(9) 



where A and B are also functions of x and Oq and include 
integrals over conformal time that trace the geodesic 
from the surface of last scattering to observation. The 
angles Or and (/>fl are not the observing angles; those are 
rather — t: ~ 9ji and (/>q|^ — tt + 4>r. The sign on the 
cos((/>ij) term (or alternatively, the (f>R to '/'obs transfor- 
mation) determines the handedness. Then a determines 
the amplitude of the fluctuation and x the pitch angle of 
the spiral. The vorticity is then 



(-) 



^2(1 + h)il + 9h) 



(10) 



Note t hat the shear and vorticity values in our original 
paper il.Taffe et alJl2f)05h contain an error in amplitude, 
although the basic conclusions are not affected. 
Equation ini can be rewritten as 



^ cx cos{(j)R ± . 



(11) 



In other words, for a given 9r, the temperature vari- 
ation follows a cos((/)i{) dependence. The phase shift </> 
is ultimately a function of 9^ and the two physical pa- 
rameters, X and Qq. The result is a spiral pattern with 
approximately N = 2/Trx twists. The smaller the x, 
the smaller the scale at which the basis vectors change 
their orientations and the tighter the resulting spiral. In 
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Fig. 1. — Examples of left-handed Bianchi anisotropy tem- 
plates in orthographic projection, all on a common color scale 
to show the relative amplitudes. These must be multiplied 
by a factor of a = {a/H)o, i.e., the shear (realistically of 
order < 10~^), in order to give the amplitude of the ob- 
served anisotropy in fiK. Note that these have been rotated 
by /3 = —90° to move the center of the structure from the —z 
pole (as defined in equation |^ to the Galactic Center. 



the case of i^o < 1 models, geodesic focusing leads to 
an asymmetry wherein the spiral structure appears com- 
pressed in one direction along the rotation axis. 



The template is calculated as 



AT 



, i.e. the contents 



(<T/mo ' 

of the curly brackets in equation times the average 
CMB temperature, so that the shear {a/H)o is the ampli- 
tude of the template to be found by fitting it against the 
CMB anisotropics, AT. Examples are shown in Figure^ 
where the template is plotted without the normalization 
by the shear. In generating all of our Bianchi templates, 
we have taken the redshift to the surface of last scatter- 
ing, or recombination, as Zrec = 1100. (Changing to, 
for example, Zrec = 1000 lowers the amplitude of the 
anisotropy by ~ 15%, implying a corresponding increase 
in the value of the shear {a/H)^ for a given AT.) 

We make the simple and pragmatic assumption that 
the anisotropy induced by the geometry simply adds to 
the statistically isotropic and Gaussian component. 

We examine a grid of such models over 0.1 < JIq < 1.0 
in increments of 0.05 and over 0.1 < a; < 10.0 in in- 
crements of 0.05 in the interval 0.1 < x < 1.0 and 
then logarithmically sampled up to a; — 10. A finer 
grid was also examined surrounding the best-fit model, 
0.52 < a; < 0.68 and 0.42 < f^o < 0.58 in increments 
of 0.02. For the largest values of x, the spiral has al- 
most disappeared (because the scale on which the ba- 
sis vectors change orientations becomes larger than the 
horizon size), and so models of higher x are self-similar. 
Smaller values of x s t art to become physically unrealistic. 
iCollins fc Haw king I \19T^ point out that for x ^ 0.05, 
the characteristic length scale over which basis vectors 
change orientation becomes comparable to the size of 
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large-scale structure, which means that lower values are 
ruled out by observations of large scale homogeneity. 
Furthermore, as discussed in i)5.3l small values of x re- 
quire higher precision analysis than is feasible. 

3.2. Data 

For this work, we arc interested only in large scale 
structure. In all of the following, unless otherwise noted, 
we use maps m HEALPix*^ format (Gorski ct al. 2005) at 
a resolution of ^gj^jg = 32 and smoothed to an effective 
beam of FWHM 5°. 5, with harmonics up to fmax = 64. 

The following full-sky maps are used in this analysis 
(where all WMAP data products are from the first-year 
data release): 

• The full-sky WMAP Internal Linear Combination 
(WILC) map release d by the WMAP team (see 
IBennett et al]l2Q03b(l . This map is formed by tak- 
ing linear combinations of the different bands such 
that the foregrounds, each of which has a different 
spectral dependence from the CMB, are removed 
leaving only CMB. The different weights of the lin- 
ear combination are determined solely by the data, 
via minimum variance, rather than by any prior 
assumptions about the foreground behavior. 

• The Lagrange Internal Linear Combination (LILC) 
map of Erikscn ct al. (2004a, 2005). The weights 
used to form the WILC map are slightly sub- 
optim al with respect t o the minimum-variance cri- 
terion l)Eriksen et al.l20 05'). and this is corrected in 
the LILC map, which uses Lagrange multipliers to 
compute the ILC weights. 

• The foreground-cleaned map of iTegmark et al.l 

(f200l), hereafter TOH. This map is also gener- 
ated by a linear combination of bands, where in 
this case, the weights are determined in harmonic 
space. 

All of these maps contain residual foreground emission, 
some of which is visible by eye along the Galactic plane 
and some of which extends to high latitudes. It should 
be noted that none of these maps is intended for high- 
precision CMB analysis, but we nevertheless use them in 
the following to locate the best-fit Bianchi template by 
full-sky convolution. Simulations show that these fits are 
affected by two opposing biases (see 5 12. 21 and S I5.3|I that 
are larger than the effects of the foreground residuals 
(see i)4.2|l . thus justifying our use of these maps despite 
their known disadvantages. In general, we use the full- 
sky maps initially to locate best-fit axis for each Bianchi 
model (see i)4.2|l . and then verify the amplitude using 
partial-sky algorithms on the following additional data: 

• WMAP uncorrected maps for each of the five fre- 
quency bands, co-added from each differencing as- 
sembly using noise weighting (see IBennett et al.l 
l2003aj) and Iso noise-weighted, coadded combina- 
tions of bands Q+Y. V+W, Q+V-l-W, Q-V, V-W, 
Q-W. 

• KpO intensity mask, excluding 23.2% of the pix- 
els in which the K-band intensity is high and also 

^ |http : //healplx ■ j pi ■ nasa . govT] 



0°.6 around known point sources, downgraded to 
^side = 32. 

Finally, we use observations at other wavelengths as fore- 
ground templates: 

• the lFinkbeiner et al.l l)1999f) model for thermal dust 
emission (hereafter FDS); 

• the iSchlegel et all l|1998t) lOO/^m intensity dust 
template (hereafter SFD), which is used a n alter - 
native to the FDS model (see discussion in S I5.2.2|I : 

• the ? Ha template, with dust correction /^j = 0.5; 

• the iDickinson et al.l (|2p03) Ha template with no 
dust correction, which is used as an alternative to 
the Finkbeiner template; 

• and the iHaslam et all l)1982() 408 MHz map of 
synch rotron emission processed by iDavies et alJ 

These foreground components are fit simultaneously to 
each band over the incomplete sky using the KpO mask, 
which reduces the effects of foreground contamination on 
the fit amplitude (see S I4.8|I . Note that although we are 
simultaneously fitting the foreground components, these 
templates are not accurate enough in the Galactic plane 
region for full-sky fits to be reliable. 

3.3. Gibbs Samples 

In addition to the WMAP data products, we also 
analyze a set of Gibbs sampled maps that were 
generated by the m ethod desc ribed bv l,Te\yell et all 
l|200irWMldelt et all lf200|) . and lEriksen et all l|2004c{r 

Effectively, this method samples the space of CMB sig- 
nal maps that are consistent with the data, taking into 
account both noise characteristics and limited sky cover- 
age. Thus, each single Gibbs sample represents a full-sky, 
noiseless CMB signal consistent with the data assuming 
Gaussianity, and the distribution of such maps describes 
the full CMB signal posterior distribution. 

Such sampled maps can thus be analyzed very effi- 
ciently using the total convolver method described above, 
since neither sky cut nor non-uniform noise complicate 
the analysis. These allow us to avoid the problem of fore- 
ground residuals in the Galactic plane, since this region 
of the Gibbs samples contains only CMB signal that is 
either consistent with the structure outside the plane, in 
the case of large enough scales, or entirely Gaussian ran- 
dom, in the case of smaller scales. The ensemble of fit 
results then reflects how well the template fits the CMB 
signal posterior distribution. In the following, we ana- 
lyze ensembles of 1000 samples corresponding to each of 
the three cosmologically important WMAP Q, V, and W 
bands. 

3.4. Assumed Signal Covariance 

Given that we are searching for evidence of anisotropy, 
the description of the expected signal covariance is not 
trivial. Bianchi models in particular are not compatible 
with inflation theory and do not make any prediction for 
fluctuations at the surface of last scattering. Clearly, a 
self-consistent theory is required to explain the observed 
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anisotropies in addition to the Bianchi component, and 
in particular, that theory must be consistent with the 
acoustic peaks now detected at smaller scales. No such 
theory currently exists, but we note that the Harrison- 
Zel'dovich power-law spectrum prediction predates infla- 
tion theory. Because it has been shown to match the 
data very well on small scales, we use the inflationary 
prediction as a starting point. 

The signal covariance expected after subtraction of any 
Bianchi component is then assumed to be that of Gaus- 
sian, isotropic CMB fluctuations fully characterized by 
the power spectrum. We use the best-fit WMAP the- 
oretical power-law spectrum to perform our fit. One 
could then refine the input spectrum based on the result 
(i.e., do a new parameter estimation using the corrected 
sky) and iterate. In the present analysis, however, we 
do not aim to improve the power spectrum estimation. 
Template fitting proves to be insensitive to the assumed 
power spectrum. (The fit result changes by less than 3% 
when using a flat, Q = 18/iK power spectrum instead.) 
So for the purposes of this analysis, the best-fit WMAP 
theoretical power-law spectrum is sufficient. 

4. PERFORMANCE, BIAS, AND ACCURACY 

In order to interpret the results of the analysis using 
real data, we need first to quantify the effects described 
above. The model selection accuracy, the bias due to the 
maximization over rotations, any bias due to foreground 
residuals, and the distribution of chance alignments are 
all effects that we can quantify using simulations. 

These are gene r ated by the LILC simulation pipeline 
of lEriksen et alJ l|2004al l2005j) . The simulations start 
with a Gaussian CMB signal generated from an as- 
sumed power spectrum and are then smoothed to the 
beam width of each WMAP differencing assembly. Pixel 
noise is added, uncorrelated and following the instru- 
ment properti e s and observation pattern described in 
iBennett et al.l l)2003a|) . Finally, the three foreground 
components above are added to create simulated raw 
data for each of the 10 differencing assemblies. The 
LILC algorithm is then used to reconstruct the corre- 
sponding processed, foreground-cleaned sky. Although 
these are known to underestimate somewhat the amount 
of residual emission along the Galactic plane, they pro- 
vide a vital indication of the morphology and approxi- 
mate amount of such residuals that may be present in 
the WILC or LILC maps. 

We apply the fitting methods outlined above to the 
ensemble of LILC simulations, with and without an ad- 
ditional known anisotropic signal, to characterize how 
well the methods perform. In most of the analysis be- 
low, a set of 1000 simulations were used in the full grid 
searches and cut-sky pixel space fitting. An expanded 
ensemble of 10,000 LILC reconstructions was used to re- 
fine the significance measures for the two best-fit models 
found as described in § |S| 

4.1. Model Selection Accuracy 

First, we add a known Bianchi component (the partic- 
ular template and amplitude found in our initial analysis 
|jaffe et al. 2005!) to a set of LILC simulations and per- 
form the full sky search over all rotations (using the total 
convolver) and over the grid of models. We find that the 
most significant model returned is close (± 0.1 in a; and 



ilo) to the correct model in ~ 50% of cases. Among the 
other ^ 50%, a qualitatively different model was found 
to be the best- fit, but the correct model was still found 
to be over 99% significant in most cases. In other words, 
only in ^ 23% of realizations was the correct model not 
detected. 

We must then see if we can distinguish the correct 
model from a false detection by other means such as in- 
complete sky fits with simultaneous foreground template 
fitting. These give an idea how much the full-sky fit is 
affected by residuals in the Galactic plane. Furthermore, 
models that appear far apart in the model space may in 
fact be fitting to the same CMB structure. We therefore 
select the several most significant models to examine in 
more detail. Then we look at what structures are driv- 
ing the fits and how they behave when the Galactic plane 
is excluded and foreground templates simultaneously fit. 
These tools give an additional qualitative way to com- 
pare different model fits. 

4.2. Full-Sky Fitting Accuracy 

Next, we consider a known Bianchi component added 
to the input noiseless, pure CMB realization (as opposed 
to the LILC reconstruction) and see how well its position 
and amplitude are recovered by the full-sky fit. For 1000 
simulations, a Bianchi component (at the same position 
and amplitude as our best-fit against the real data) is 
added to the input CMB sky and then fit using the total 
convolver method described above. In ^ 80% of realiza- 
tions, the returned fit is within 5° (approximately the 
beam width) of the correct location. (In the orientation 
angle, it is less accurate due to the self-similarity of the 
spiral structure under such rotations. The returned ori- 
entation is within 10° in 52% of the simulations.) The 
amplitudes average ~ 7% higher than the input value (as 
noted in i l2.2(l . with an rms error of about 80% the cal- 
culated error. Neither of these facts is unexpected, since 
these values are the selected maxima, and their distribu- 
tion is not Gaussian. The results are quantitatively the 
same for the LILC reconstructed skies, indicating that 
the foreground residuals do not introduce a significant 
additional bias in the case in which a real Bianchi com- 
ponent is being fitted. Note that simulations in which 
the input Bianchi model has an amplitude a factor of 
^ 3 higher show a much smaller relative bias (~ 1%), as 
one would expect. 

4.3. Cut-Sky Fitting Accuracy 

The cut-sky fits are performed with the Bianchi model 
at the fixed location found as the best-fit using the full- 
sky total convolver method. As described in § 12.21 there 
is a bias introduced by the selection of the maximum am- 
plitude position. This bias will also be reflected in the 
cut-sky fits, although masking out the Galactic plane 
should remove some of the bias due to residual fore- 
ground emission. 

For fits to the raw data outside the KpO mask, eight 
template components are fit simultaneously to each 
band, the three foreground templates described in S I3.2I a 
monopole term, the three spherical harmonics represent- 
ing the real- valued dipole terms, and the Bianchi tem- 
plate. 

For simulations with no additional Bianchi component, 
the results show amplitudes on average 6% lower than 
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Fig. 2. — Distributions of fit results for the Bianchi compo- 
nent for 1000 simulations without {black lines) and with (red 
lines) a Bianchi component added. Vertical solid lines show 
the means, and vertical dashed show the actual errors. The 
vertical green line shows the true value and expected errors. 



those from the LILC fits. This is further indication that 
chance alignments are affected by residuals in the plane, 
since the exclusion of that region tends to lower the fit 
amplitude. 

Simulations with an additional Bianchi component at 
a known position and amplitude were run through the 
same pipeline, i.e., first the full-sky LILC reconstruction 
was used to find the best-fit location, then that location 
used to fit the template to the cut sky in pixel space. As 
described above, the total convolver will return a position 
that is very close to the true position but one where the 
fit amplitude happens to be highest due to CMB and 
noise contributions. These will also affect the cut-sky 
fits, which also show a bias of 3%. This is lower than 
the bias in the full-sky fits, showing that a few percent of 
the full-sky bias is due to residuals in the Galactic plane 
region. The relative drop in amplitude between the fuU- 
and cut-sky fits for true detections is on average half the 
drop in the case of chance alignment detections. 

Figure |3 shows what these distributions look like for 
the fit to 1000 simulations in the V band, both in the case 
where a Bianchi component is added {red histogram) and 
where it is not {black histogram). Also plotted as vertical 
lines are the mean and rms errors on the distributions, 
and the true value and expected errors plotted in green. 
The small bias in the value of the Bianchi fit is seen in 
the distance between the vertical red and green lines. 

Note that in all these cases, the bias in the fits affects 
the absolute amplitude (i.e., shear) estimate, but not the 
significance of the fit, since the ensemble of simulations 
used to estimate the significance is also affected by such 
a bias. The expected bias in the amplitude is also much 
smaller than the error bar. Therefore this does not affect 
our basic results, namely the particular best-fit model, its 
location, its approximate amplitude, and its approximate 
significance relative to chance alignments. 

4.4. Chance Alignments 

For a given sky realization, we find the best model as 
described in § 12.31 and then simply compare the ampli- 
tude of that fit against the ensemble of amplitudes for 
that model relative to Gaussian simulations to estimate 



the significance. Visual inspection of the WMAP sky 
maps shows no obvious Bianchi component, so any such 
signal must remain at or below the level of the stochas- 
tic component. Chance alignments may therefore either 
cancel a Bianchi-induced signal or give a false positive . 
The former effect was quantified in S 14. II at ~ 23%, but 
the latter is more difficult to quantify. 

The family- wise error rate (FWER), the expected 
number of false detections when testing m hypotheses, 
is ^ mp when p is the probability of one false detection. 
If 3(7 is the detection threshold (implying p ~ 0.003) 
and one tests 100 different hypotheses (or models), the 
FWER is then 0.3, meaning one gets a false detection 
somewhere in the model space one-third of the time. 
Over our grid of Bianchi parameters, the models are not 
independent (since models close in {x,ilo) space will re- 
semble each other closely) , so we cannot determine a pri- 
ori what the true frequency of false detections would be, 
but we can get this from the ensemble of Gaussian sim- 
ulations. 

We perform the full-sky search using the total con- 
volver over the grid of Bianchi models and find the best- 
fit model for each realization. We find that a false detec- 
tion due to a chance alignment that has a significance of 
99.7% occurs in ^ 17% of the cases. A better compari- 
son might be to use the representing the goodness of 
the fit. We then co mpa re the statistic rji = {xq — Xi)/Xo 
(defined above in ii2.3|l , namely the relative improve- 
ment in the when the Bianchi model is subtracted. 
We find that by this measure ~ 10% of the best chance 
alignments fit their respective realizations as well as our 
best-fit model does the WMAP data (see §0). Note, 
however, that these statistics are dependent on the as- 
sumed amount of large-scale power. The above numbers 
simply imply that a detection of a Bianchi model with 
an amplitude higher than in 99.7% of simulations is more 
than 4 times as likely to be real as it is to be a chance 
alignment, in the absence of all other information. 

5. APPLICATION TO THE FIRST- YEAR WMAP DATA 

Armed with the information gained from the analyses 
of simulations, we can now examine the fits to the real 
data. 

5.1. Fits Over Model Space Grid 

Using the total convolver to find the best orientation, 
we fit the grid of Bianchi models to each of the WILC, 
LILC, and TOH full-sky processed maps. Figure El shows 
filled contours over this grid for the LILC. (The results 
for the WILC and TOH look very similar.) For each 
point on the grid corresponding to a model of the given 
{x,Qo)i the template is fit to the LILC map, and the 
color indicates the significance estimate of the resulting 
amplitude, i.e., the fractional number of Gaussian LILC 
simulations (out of 1000) with lower amplitude. As dis- 
cussed in § 12.31 we use a finer grid and better method 
for determining the best-fit model and there by select a 
slightl y different model than the analysis in l.laffe et alJ 
l)2005ft . But it is apparent from the right panels of Fig- 
ure O that the significance as a function of the Bianchi 
parameters x and flo is fiat in the region ±0.1 in both x 
and f2o about the maximum. 

We find that the most significant fit is found with 
a right-handed Bianchi template of a: = 0.62 and 
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Fig. 3. — Significance as percentage of LILC simulations 
whose best-fit cliance alignment amplitude is lower. Left 
panels show left-handed models, and right panels show right- 
handed models. Over plotted contours are at 99.3%, 99.5%, 
99.7%, and 99.9%. Two color scales are used to show the 
global structure ([emphatop panels) as well as that near the 
peaks [bottom panels). 



r^o = 0.5 when that template is rotated to a posi- 
tion and orientation given by Euler angles (following 
the total convolver's "zyz" convention about fixed axes) 
($2,6,*i) = (42°, 28°, -51°). As defined in § O the 
spiral structure of the unrotated model is centered on the 
south pole (or —z axis), so this rotation places the cen- 
ter of that structure at Galactic longitude and latitude of 
(Z,6) = (222°, —62°) and changes it's orientation about 
that location by <i>i — —51°. This model fits at an ampli- 
tude of (f-)g = 4.29 X 10~i^ which is higher than 99.7% 
of the 10 000 simulat ions. This mode l and the best-fit 
from previous work llJaffe et alJl27iol at a; = 0.55 are 
almost identical. 

All models near this best-fit (x, fig) return the same 
location for the center of the spiral within 3° but vary 
the orientation (Euler angle $1) up to 36°. The broad 
spiral in all of these models is very self-similar under 
these rotations, so the change is driven largely by the 
precise locations of the paired hot and cold spots. 

Looking at Figure |31 one can see that more than one 
model appears "significant" in the sense of fitting with 
an amplitude above 99% of the amplitudes found fitting 
that same model to Gaussian simulations. As discussed 
above in t l4.1l this is not surprising, and we must examine 
each of these models in more detail. 

The full resolution ILC map is shown along with the 
best-fit Bianchi model on the same scale and the cor- 
rected ILC map in Figure^ A summary of all fit results 
is shown in Table The expected bias in these results 
is discussed in ii4.3l The following sections describe the 
two most interesting models in more detail. 

5.2. Two Best Fits 

5.2.1. Left-handed Model {x^Uq) = (0.62,0.15) 

The most significant left-handed model, at 99.4%, is 
at (x, fio) — (0.6 2,0.15). This m odel was not found in 
our earlier work l)Jaffe et al.ll2005j ). because it is only in 



TABLE 1 

Fitted template amplitudes 







{w/H)o ^'dosiml < l"obsl) 


Map 


(xlO"!") 


(xio-i") % 



Right-handed (a;,Qo) = (0.62,0.5) 



WILC 


4.33 ±0.82 


9.58 


99.8 


LILC 


4.29 ±0.82 


9.49 


99.7 


TOH 


4.03 ±0.82 


8.92 


98.6 




2.59(4.13) ±0.83 


5.72 


16.7(99.1) 


Ka'' 


3.50(4.09) ±0.83 


7.74 


86.9(99.0) 




3. 76(4. 11) ±0.83 


8.31 


95.6(99.1) 


ya 


3.99(4.19) ±0.83 


8.82 


98.1(99.5) 




4.08(4.35) ±0.82 


9.03 


99.1(99.8) 


QVW 


3.84(4.15) ±0.83 


8.49 


96.8(99.2) 


VW=' 


3.99(4.22) ±0.83 


8.84 


98.2(99.6) 




0.06(0.11) ±0.02 


0.13 


99.0(100.0) 




-0.05(-0.08) ± 0.02 


0.11 


93.8(99.0) 


Q-W 


0.01(0.04) ±0.02 


0.02 


25.0(83.2) 




4.09 ±0.10'= 


9.04 






4.11 ±0.10'= 


9.09 






4.12 ±0.11'= 


9.12 





Left-handed (a', Ho) = (0.62,0.15) 



Note. — Amplitudes of the best-fit model derived from various 
combinations of data and various methods as described in the 
text. The full sky was used in the analysis of the WILC, LILC, 
TOH, and Gibbs samples, while the KpO mask was imposed for 
the remaining maps. 

^Simultaneous fits with foreground components. In parentheses 
arc the values using the SFD dust template instead of the FDS, 
and the Dickinson ct al. 1 20^) Ha instead of ?. 

''Average over 1000 Gibbs samples. 

'^Errors are rms variation over Gibbs samples. 



a fairly small region of the model space that this fits 
with any significance, and our previous, coarser grid ef- 
fectively straddled the peak in The best- fit loca- 
tion for this model puts the center of the structure at 
(/, b) — (320, —20), which is closer to the Galactic Center 
region than the best-fit right-handed model, raising the 
question of how much it is driven by foreground residuals. 

Cut-sky fits give fit amplitudes for this component that 
are 8% lower and significances of ~ 96% in most cases. 
Furthermore, the Galactic center region tends to draw 
the template in simulations; the best-fit location among 
the simulated LILC maps for this model is twice as likely 
to be found in the area around (0°, —20°) as should be 
expected from a uniform distribution. The only thing 



WILC 


2.39 ±0.47 


22.31 


97.8 


LILC 


2.49 ±0.47 


23.29 


99.4 


TOH 


2.45 ± 0.47 


22.94 


99.0 




2.33(3.31) ±0.50 


21.76 


96.3(99.9) 


Ka^' 


2.24(2.63) ±0.50 


20.93 


94.8(99.3) 




2.29(2.50) ±0.50 


21.42 


96.0(98.9) 


ya 


2.33(2.44) ±0.50 


21.81 


96.7(98.6) 




2.32(2.46) ±0.49 


21.69 


96.3(98.4) 


QVW^ 


2.30(2.48) ±0.50 


21.46 


96.1(98.8) 


VW» 


2.34(2.44) ±0.50 


21.85 


96.7(98.6) 




0.03(0.02) ±0.02 


0.27 


78.6(63.1) 


V-W'' 


-0.06(-0.10) ± 0.02 


0.55 


96.4(99.9) 


Q-W^ 


-0.03(-0.08) ± 0.02 


0.29 


73.1(99.3) 


Q" 


2.10 ±0.11= 


19.67 






2.08 ±0.11= 


19.46 




W*" 


2.09 ± 0.09'= 


19.53 
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Fig. 4. — Top: WMAP Internal Linear Combination map. 
Middle: Best-fit Bianchi Vllh template (enhanced by a factor 
of 4 to bring out structure). Bottom: Difference between 
WILC and best-fit Bianchi template; the "Bianchi-corrected" 
ILC map. Over-plotted on each as a dotted line is the equator 
in the reference frame that maximizes the power asymmetry 
as described in § 16.31 



that all of the LILC simulations have in common is fore- 
grounds, so this is an indication that there is some resid- 
ual there that is a weak attractor. One possibility is 
the "free-free haze" described by Finkbeiner (2004, see 
also Patanchon et al. 2005), although this haze does not 
match up well with the template structure, the two show 
little cross-correlation, and inclusion of Finkbeiner's haze 
template in the simultaneous fitting does not alter the fit 
amplitude of the Bianchi model. 

In FigureEl it looks like the fit should be largely driven 
by the cold region below the Galactic center. The cross- 
correlation maps described in § 12.41 do show correlation 
there but also indicate that the fit is largely driven by a 
very strong signal in the Galactic plane. Figure 151 shows 
these maps for both this model and the best-fit right 
handed model. Where the right handed model shows 
relatively uniform correlation over the hemisphere about 
the best-fit axis, this model shows a rather concentrated 
region including a very strong driver on the Galactic 
plane. 

The Gibbs samples throw further doubt on this model. 
Among the 1000 Gibbs samples in each of Q, V, and W 
bands, this model fits at the same approximate location 
as for the LILC map less than half of the time. Where 
the location was the same, the amplitude of the best-fit 
is significantly lower for the ensemble of Gibbs-sampled 
maps, which drop over 15% in amplitude to a mean of 
2.1 X 10^^", indicating that some of the structure in the 
data that drives the fits is not consistent with the poste- 
rior CMB distribution as determined by the Gibbs sam- 
pling technique. Furthermore, this model is almost as 
likely to fit near the location of the best-fit right-handed 
model instead of near the Galactic center. This is largely 




Fig. 5. — Two significant models. The left panels show the 
best-fit left handed model with {x,Q,o) = (0.62,0.15), while 
the right panels show the best-fit overall model, right-handed 
with {x,Q.q) = (0.62,0.5). The top panels show the template 
amplified by a factor of three to bring out the structure. The 
middle panels show the corresponding cross correlation map 
(see H2.4II scaled from —1% to 2%. The bottom panels show 
the "corrected" WMAP Q-l-V-l-W map scaled from -150 to 
150/i-ft'. The grey region is the excluded region of the KpO 
mask. 



driven by the cold spot. 

In summary, this model is quantitatively less signifi- 
cant than the best-fit right-handed model based on the 
cut sky and Gibbs sample fit values. Furthermore, the 
morphology indicates that foreground residuals drive the 
full-sky fit. 

5.2.2. Right-handed Model {x,Q.o) = (0.62,0.5) 

Figure 13 shows that the best-fit model is this right- 
handed model. 

The amplitude of the best-fit Bianchi component varies 
somewhat across the different frequencies, in all cases 
lower than the full-sky amplitude fit with the LILC. As 
discussed in t |4.3l this is likely due to small foreground 
residuals, but does not mean that the detection is a false 
positive; the same effect occurs in simulations that in- 
clude a Bianchi component. The amplitude in the W 
band, in which the least foreground residuals are ex- 
pected, is still higher than ^ 99% of simulations. The K 
and Ka band fits are significantly lower when the FDS 
dust and Finkbeiner Ha templates are used, but are con- 
sistent with the other bands when the SFD dust and 
Dickinson Ha templates are used instead. It is known 
that foreground subtraction is a problem even at high 
latitudes in the K and Ka bands, and these residuals are 
clearly affecting the low frequency fits. Looking at the 
residuals of the two fits shows that the difference may 
be driven by a small region around {I, b) = (300°, —15°) 
where the dust templates differ strongly. The higher fre- 
quency fits, however, are more consistent. The difference 
maps, e.g., Q-V, should contain no CMB component but 
only foreground residuals and noise. The fact that the 
Bianchi component amplitude found from these maps is 
less than 2% of the co-added map amplitude is an indica- 
tion that such residuals are not contributing significantly 
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to the fit. 

The results of fitting the Gibbs-sampled maps show 
that for this model, the amphtude is quite stable over 
the ensemble of Gibbs samples, with, e.g., a mean of 
(4.12 ± 0.1) X 10-1° in the W band compared to 4.08 x 
10-1° for the cut-sky fit to the raw data. Since the Gibbs 
samples represent the posterior CMB distribution, taking 
into account foregrounds and iterating over the power 
spectrum, these results are a strong indication that the 
fit is due primarily to CMB signal. 

Figure [S] shows the cross-correlation map as described 
in § 12.41 which give a visual indication of what regions 
drive the fit. Unlike the left-handed model (left), which 
shows one concentrated region in the Galactic plane to 
be driving the fit, this model correlates over more than 
half the sky at moderate levels. One can see that the 
cold spot does partly drive the fit, but no particular 
region can be said to dominate. Fits to the combined 
QVW and VW data where the cold spot is excluded (in 
a 10° radius around {l,b) — (209,-57)) have compara- 
ble amplitudes to fits where the region is included (only 
6% lower) . Further Gibbs samples were also computed 
while masking this region. Full sky searches using these 
samples show that fewer than 20% return positions more 
than 10° from the original location, and amplitudes that 
are on average 15% lower (which is within the calculated 
error bar). These results confirm that the cold spot does 
affect but does not exclusively drive the fit amplitude. 

5.3. Location and Orientation Accuracy 

As mentioned above, where a Bianchi component was 
added to simulations at a known location, the full-sky 
search with the total convolver returned the correct posi- 
tion within 5° in ^ 80% of realizations. The uncertainty 
in the location is due to the CMB fluctuations, which 
are quite comparable to the Bianchi component at the 
amplitude detected. 

To determine how the amplitude changes with the po- 
sition and orientation of the template compared to the 
data, we take the best-fit Bianchi model and fit it to the 
LILC on a grid of fixed positions within 20° of the best-fit 
position. Results are shown in Figures IHl and The ori- 
entation is not very sensitive in this model, whose spiral 
structure is self similar under rotations about its symme- 
try axis; only the precise positions of the hot and cold 
spots affect the variation with orientation angle. The am- 
plitude drops by 1% when the orientation is 4° off. The 
location of the symmetry axis is a bit more sensitive, 
where the amplitude drops by 3% at 2°. The fact that 
the total convolver at this resolution uses steps of 2°. 8 
means that its best-fit amplitude can be several percent 
off of the actual maximum. All the fits to the simulations 
as well as the data are subject to this same uncertainty. 
If we assume the worst, that the LILC amplitude was 
found at its true maximum (i.e., the true axis of sym- 
metry happened to lie exactly on the center of one of 
the total convolver's bins) and the simulations are all 
at 1°.4 away from their true maxima (i.e., the axis ex- 
actly between bins) and have true values correspondingly 
higher, the comparative significance could then be over- 
estimated by only 0.5%. The likely effect is of course 
much smaller. 

Using the LILC map at higher resolution, £max = 128, 
gives an accuracy in the total convolver of 7r/£niax = 



distance (degrees) 



Fig. 6. — Average fit amplitude as the location of the tem- 
plate varies from the best-fit position. For these fits, the 
orientation of the template is unchanged. 



Euler orientcition angle (degreeK) 

Fig. 7. — Fit amplitude as the orientation of the template, 
i.e., the Euler angle 7, is changed. For these fits, the location 
in longitude and latitude is unchanged. Note that this grid 
finds a preferred orientation 2° from that found by the total 
convolver (due to the slightly different grids used.) 



1°.4. The position returned is identical, with only the 
orientation one step of 1°.4 different. 

The above apphes to the best-fit model at (a:, rig) — 
(0.62,0.5), but other models have structure at different 
angular scales. In particular, for the region of small x and 
r^o, where a tightly wound spiral is even more tightly fo- 
cused in one hemisphere, the fit amplitudes are far more 
dependent on the exact position. Because the total con- 
volver resolution is 7r/£niaxj our analysis is not as sensi- 
tive for this region of model space as it would be for a 
higher resolution analysis. In these cases, the difference 
of a few degrees can mean a large difference in amplitude. 
Simulations show that, although the location returned is 
the closest bin to the true location, the amplitude of a 
model {x, fto) = (0.1, 0.1) is underestimated by ~ 20% on 
average due to the limited resolution. Increasing the res- 
olution of the analysis to HEALPix A^gj^jg — 64 increases 
the mean and brings it closer to the correct value, but it 
is still underestimated. (Higher resolution analysis with 
the total convolver is not feasible due to the memory and 
CPU requirements.) In the region of model space where 
X > 0.25 and fio > 0.25, this effect drops to less than a 
few percent. 

A more detailed look at these models at increased res- 
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olution (-^side ~ shows no evidence that the lower 
resolution analysis missed a significant detection. But 
the limits placed on shear and rotation are less stringent 
than they would be were a higher resolution analysis fea- 
sible. 

5.4. DMR Fit 

Our best-fit amplitude is below the upper limit DMR 
could place on the shear. Using this model, a fit to the 
DMR data gives {jj)^ = 3.38 ± .98 x 10-^°, which is 
within our best-fit error bar for the WMAP data, but 
which is not distinguishab le from a chance alignment for 
DMR. lKogut et all l)1997|) report a distribution of T val- 
ues for chance alignments up to 4.5. Our fit value and 
error give T — 3.4, and although this value comes from 
different methods and assumptions, it is roughly compa- 
rable. 

5.5. Sensitivity to Assumed Power Spectrum 

As mentioned in § 13.41 assumptions about the cosmo- 
logical parameters go into this analysis from the begin- 
ning with the choice of the signal covariance matrix. In 
effect, we are assuming that the CMB signal consists of 
an anisotropic Bianchi-induced component plus a statis- 
tically isotropic, Gaussian random field described com- 
pletely by its power spectrum, which is taken to be the 
WMAP best-fit theoretical power law spectrum. As we 
are searching for evidence of a model that affects the 
power spectrum at large scales and that is inconsistent 
with infiation, this approach obviously lacks consistency. 

We have verified, however, that changing the assumed 
parameters and using, for example, a flat (5 = 18 mK 
power law spectrum, or a completely implausible spec- 
trum, has little effect (less than 3%) on the resulting 
best-fit amplitude and position for the Bianchi compo- 
nent. In fact, the power spectrum affects only the es- 
timated significance of the result, as that significance 
is dependent on the expected level of large-scale CMB 
structure that drives chance alignments. As shown in 
Figure |H1 correction for this Bianchi model lowers the 
large-scale power. Our significance estimates are based 
on simulations generated assuming a higher level of large 
scale power, so the significance of the detection would in- 
crease when compared to an ensemble consistent with the 
corrected power spectrum. 

6. IMPLICATIONS 

There are several interesting results based on WMAP 
first-year data that are inconsistent with the assump- 
tions of isotropy and Gaussianity and that are imme- 
diately relevant to this st udy. De Ol iveira-C osta ct al. 
f2nn4V lLaud fc Ma.gueiinl (j^^), a,nd lCnni et a,] J p.m(h 
(and sources therein) examine the low-£ multipoles of 
the foreground-cleaned data and find that, in addition to 
the anomalously low quadrupole amplitude, the preferred 
axes of the quadrupole and octopole are anomalously 
weU aligned in the d irecti on of (I, b) = (-110° 60°). 
lEriksen et aP l|20n4bf) and iHansen et all l|2n04br> find 
a system of reference (roughly aligned with the eclip- 
tic) in which there is a significant difference in large- 
scale power between the two hemispheres at the 98%- 
99% level, w ith significa ntly more p ower in the south. 
iVielva et all p0041 and iCruz et all (,2005) detect non- 
Gaussianity in the WMAP combined Q-V-W map using 
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Fig. 8. — Comparison of power spectra. The gray and black 
solid lines show the power spectrum estimated from the co- 
added V-l-W map before and after correcting for the Bianchi 
template, respectively. The dotted gray and black lines shows 
the theoreti cal best-fit power-spec tra from the WMAP-team 
analysis and IHansen et ahl i2004afl respectively. The latter is 
a fit to the northern hemisphere data alone. The dashed grey 
line is the power in the Bianchi template alone. 



spherical wavelets; they find significant kurtosis in the 
wavelet coefficients at a scale of 10° and identify a cold 
spot at (Z, b) = (209°, —57°) as the probable source. Our 
choice of models was partly motivated by the morphol- 
ogy of these anomalies, and indeed, subtracting for our 
best-fit Bianchi template corrects them. 

6.1. Quadrupole Amplitude 

The quadrupole amplitud e has been considered anoma- 
lously low since COBE Csee lde Oliveira-Costa et al.l2004] 
and references therein). As pointed out bv IJaffe eT'aLl 
(2005), the correction for this Bianchi component raises 
the low quadrupole amplitude to a value more consistent 
with the theoretical power spectrum. This result is un- 
changed with the best-fit model of this work, since the 
models are almost the same. Should this be considered 
"fine tuning"? 

We can simulate the situation by taking as the the 
primordial q uadrupole the W MAP quadrupole (as de- 
rived bv IBielewicz et al.ll200^ minus the quadrupole of 
our best-fit Bianchi model. If we then add the Bianchi 
quadrupole at random orientations, we can see how likely 
it is that the resulting total quadrupole be as low as the 
observed WMAP quadrupole. We find that the likeli- 
hood is ~ 5%. This implies that the level of "fine tuning" 
required to end up with the low observed quadrupole is 
not exceptional. 

We further take a set of 1000 simulated Gaussian CMB 
skies, with and without a Bianchi component, and fit 
our best-fit Bianchi model to them. The "corrected" 
quadrupole is on average 5% lower than the original, 
which is to be expected considering that the fit is a least 
squares solution. In contrast, using the real LILC data, 
the correction has the effect of raising the quadrupole. 
This happens in over ~ 20% of the simulations, so while 
this is not the average behavior, it is not extraordinary. 

6.2. Low-i Alignment and Planarity 

De Olivei ra-Costa et al. ( 2004), iLand fc Magueiiol 
()2005(1 . and iCopi et all ()2006ri discuss the statistically 
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Fig. 9. — Low-£ multipoles of the WILC corrected {bottom 
panels) and uncorrected {top panels) for the Bianchi compo- 
nent. 



anomalous alignment of the quadrupole and octopole in 
the WMAP data. The preferred axes of the i = 2 and 
^ = 3 multipoles are only 7° apart (roughly in the direc- 
tion of {l,b) — (—110°, 60°)), which is anomalous at the 
99.3% level compared to simulations. After subtracting 
the best-fit Bianchi template, these axes lie 74° apart, 
consistent (at 27%) with the statistically isotropic simu- 
lations (see Fig. ISJ. 

The planarity of the low-.^ multipoles has 
also been considered somewhat an omalous (see 
de Oliveira-Costa et alJ 120041 and iLand & M agueiid 
2005 for a discu ssion) . Th e i-statistic defined by 
ie Olivcira-Costa et _al.l 12Q0J1 ) provides a measure of 
this planarity. Again, subtracting the Bianchi template 
lowers the significance of the the low-^ multipoles. The 
planarity of the octopole in particular drops from a 
significance of ~ 90% (depending on whether the WILC 
or LILC is used) to ~ 50%. FigureslHl (b) and (d) shows 
how the planarity of the octopole is disrupted. This will 
also impa ct the result s of multipole vector analyses such 
as that of iCopi et all ((2006). 

6.3. Large-Scale Power Asymmetry 

lEriksen et all l|2004b(l and iHansen et all l)2004bfl re- 
ported that the large-scale power {£ < 40) in the WMAP 
data is anisotropically distributed over two opposing 
hemispheres (in the reference frame in which the z-axis 
points toward {l,b) = (57°, 10°); see Fig. Unil, with a 
significance of 3cr compared with simulations. Repeat- 
ing the analysis and adopting the Kp2 sky coverage, we 
compare the corrected V+W WMAP map with 2048 
simulations. We find that ^ 14% of the simulations 
have a larger maximum power asymmetry ratio than the 
Bianchi-corrected map, whereas only 0.7% have a larger 
ratio than the uncorrected data (see Fig. I10() . It is ap- 
parent that the maximum power ratio between any two 
hemispheres is significantly suppressed after subtracting 
the Bianchi template, as no asymmetry axis is found at 
any statistically significant level. It is apparent from that 
figure, however, that some residual power asymmetry re- 
mains. This comes largely from the range 20 < I < 40, 
where the Bianchi template has little power, indicating 
that a model with more small-scale structure may be 
needed. 

6.4. Wavelet Kurtosis 
iVielva et all l|2004fl and iCruz et all l|2005(l used a 




Fig. 10. — Power ratio between hemispheres in WMAP ILC, 
corrected {bottom) and uncorrected {bottom) for the best-fit 
Bianchi component. 




Wavelet scale (deg) 

Fig. 11. — Kurtosis in wavelet coefficients. The boxes and 
crosses show the kurtosis before and after subtracting the 
Bianchi template, respectively, computed from the southern 
{dotted line) and northern {solid line) Galactic hemispheres. 



wavelet technique to detect an excess of kurtosis in the 
wavelet coefficients and isolate an unusually cold spot 
{^ 3ct significance relative to Gaussian simulations) at 
Galactic coordinates {l,b) = (209°, -57°). Referring 
again to Figure ^ we see that a cold spot is indeed 
present at the right location, in the form of the center of 
the spiral. ^^^^^^^^ 
We therefore also repeat the analysis of IVielva et alJ 
|2004), and compute the kurtosis of the wavelet coeffi- 
cients as a function of scale from both the WILC and 
the corresponding Bianchi-subtracted map. A |6| < 20° 
galactic cut is imposed in this case, for computational 
convenience. 

The results from this exercise are reported in Figure lTTI 
After subtracting the Bianchi template, the significance 
of the southern hemisphere anomaly is greatly reduced, 
and no new non-Gaussian features have been introduced. 

7. DISCUSSION AND CONCLUSIONS 
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We have considered a fast and efficient method for fit- 
ting a template to the fuU sky in harmonic space and find- 
ing the best-fit location and orientation. The total con- 
volver algorithm evaluates the correlation between the 
sky and the template at every possible relative orienta- 
tion using fast Fourier transforms. With this algorithm, 
the search for the best-fit becomes 2 orders of magnitude 
faster than the corresponding search performed one ro- 
tation at a time. This method, along with pixel-space 
simultaneous foreground fitting, provides a powerful tool 
for testing any deterministic model for anisotropy in the 
CMB. Simulations generated by the LILC pipeline allow 
us to quantify the bias, investigate the effects of fore- 
ground contaminants, and show how well each of these 
methods detects a known input. 

We have applied this method to the first-year WMAP 
data to search for evidence of shear and vorticity us- 
ing templates derived for Bianchi type VII/i universes. 
We find a surprisingly significant correlation between the 
WMAP data and a right-handed Bianchi model with 
X = 0.62, f7o = 0.5, and shear of (f)^ = 4.3±0.8 x IQ-i", 
implying a vorticity of (■^)q — 9.5 x 10^^°. The cen- 
ter of the spiral structure lies at approximately (Z, &) = 
(222°,— 62°). Simulations show that this amplitude is 
likely to be biased by ~ 7%, implying a true amplitude 
closer to 4.0 x 10~^°. Incomplete sky fits, simultaneous 
foreground fitting, and fits to a set of Gibbs samples are 
all consistent with this amplitude and indicate that con- 
fusion with Galactic emission is unlikely to contribute 
significantly to this detection. 

Correcting the WMAP data for the effect of the best- 
fit model solves several problems seen in the data. The 
corrected maps show significantly reduced power asym- 
metry between any two hemispheres. The correction also 
eliminates the non- Gaussian kurtosis in the wavelet co- 
efficients detected bv lVielva et all l)2004|) and lCruz et al.l 
((1)05), raises the low measured quadrupole by a factor 
of 2, and disrupts the planarity of the octopole and its 
anomalous alignment with the quadrupole. In short, the 
data appear far more Gaussian and isotropic after cor- 
rection. 

The orig inal analyses by iKogut et alJ l)1997t) and 
IBunn et all llT99l were limited by the signal-to-noise 
ratio level in the DMR instrument. Our best-fit result 
is just under their upper limit but still significant due 
to WMAP's greatly improved signal-to-noise ratio. Fur- 
thermore, the Kogut analysis searched a coarse (~ 10°) 
grid of possible locations and orientations, while with the 
total convolver, we can efficiently search a finer grid. 

How likely is it that our best-fit model is a true detec- 
tion rather than a chance alignment? Considering the 
best-fit model by itself and comparing its fit amplitude 
to simulations, it is higher than 99.7% of simulations. 
However, the simulations also show that 10% — 20% of 
Gaussian, statistically isotropic skies will have one of 
the Bianchi models appear as significant. Considering 
the fact that the sky is approximately Gaussian and 
isotropic, one would not expect to find a more definitive 
detection based on template fitting alone. But the dis- 
tribution of chance alignments in the simulations is sen- 
sitive to the amount of large-scale power assumed, and 
that is lowered by the Bianchi correction to the WMAP 
data. Furthermore, the cumulative probability that a 



chance alignment not only fits at the level of our best- 
fit model but also has the effect of resolving the several 
anomalies in the data must also be considered in any 
qualitative judgment of the significance of this result. 

Further improvement to the data will not refine these 
measures significantly, because at the WMAP sensitiv- 
ity level, the analysis is already very close to the ex- 
pected distribution of chance alignments in the absence 
of noise. Improved foreground subtraction will, how- 
ever, remove some of the possible confusion and bias, 
but neither higher resolution nor higher signal-to-noise 
ratio data should change this result nor be able to pro- 
vide additional information concerning the question of 
whether the fit is a real detection of vorticity and shear. 
Answering that question will require additional verifiable 
predictions for the effects of vorticity and shear on other 
observables. 

However, in the context of the anomalies that this hy- 
pothesis can explain, the possible detection is certainly 
provocative. The most important result of this analy- 
sis is that a model with vorticity and shear can explain 
the observed asymmetry in the CMB anisotropies and 
the non-Gaussian cold spot. Note that this asymmetry 
exists only in the fio < 1 versions of these Bianchi mod- 
els. Significant evidence currently indicates that is 
very close to 1, so our best-fit model cannot be consid- 
ered physically realistic. However, as mentioned in S I3.1I 
Barrow et al. (1985) did not include any dark energy 
component. Furthermore, the Bianchi model does not 
include a mechanism to generate structure at the surface 
of last scattering. A self-consistent theory is required 
that can explain the small scale fiuctuations, and in par- 
ticular the acoustic peaks, in the context of a Universe 
with shear and vorticity. But from a pragmatic point of 
view, one can conclude that, regardless of the viability of 
the particular Bianchi model, this result gives a measure 
of the significant deviation from isotropy in the data. 

We consider this result to be further motivation for 
considering ideas outside of the so-called concordance 
model of cosmology. There are anomalies in the data 
that are inconsistent with the theory of a Gaussian, sta- 
tistically isotropic universe, and Bianchi models are only 
one such anisotropic model that merits investigation. We 
have demonstrated a method of template fitting that can 
be applied to test any model that makes a deterministic 
prediction for an anisotropy pattern in the CMB. The 
best-fit Bianchi model provides a template temperature 
pattern that can explain the observed anomalies in the 
data and that describes the morphology theorists may 
need to reproduce in considering alternatives to the stan- 
dard cosmological model. 
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